#include "uhm.hxx"

#define N 1

using namespace uhm;

typedef ID_Node_H_<N>          Id_Nd;
typedef ID_Element_<N>         Id_El;

//typedef LU_nopiv_Linal_        Uhm;
typedef LU_piv_Linal_        Uhm;
//typedef LU_incpiv_Linal_        Uhm;

typedef Direct_Solver_<Id_Nd,Id_El,Uhm> Solver;

int main (int argc, char **argv) {

  initialize(argc, argv);

  // ** Environments
  Int_
    type = Mat_Base_::DOUBLE,
    uplo = Mat_Base_::FULL,
    sym  = Mat_Base_::UNSYM,
    rhs  = 1,
    blk  = 256;

  Int_
    kind = Mesh_Base_::NODE_DEFAULT;

  // ** Solver initialization
  Solver solver(type, uplo, sym, rhs, blk);

  // ** Import a mesh from a file_in
  option_begin();
  std::string fin;
  get_option_string("-in", "Reading a file", __FILE__,
                    "../../../uhmfile/toy_1.uhm", fin);
  Int_ time;
  get_option_int("-time", "Time tracking level 1) User 2) Front End 3) Internal",
                 __FILE__,UHM_TIME_LEVEL_FRONT,time);
  option_end();

  time_begin(time);

  std::cout << "Opening file ... " << fin << std::endl;
  std::ifstream in;
  in.open(fin.c_str());
  stream2mesh<Solver::_Mesh_>(in, solver.mesh());
  in.close();
  std::cout << "Closing file ... " << fin << std::endl;

  // ** Create random matrices in work
  UHM_ERROR((mesh2cache<Solver::_Mesh_,Solver::_Cache_>(solver.mesh(), solver.work())),
            ">> Error in mesh2cache");
  
  solver.work().set_matrices(type, rhs, blk);
  solver.work().create_matrices();
  solver.work().random_matrices();

  // ** create tree and allocation
  UHM_ERROR(solver.create(),    ">> Error in Creating");

  time_in("direct:: Decomposing");
  UHM_ERROR(solver.decompose(), ">> Error in Decomposing");
  Double_ t = time_out();

  solver.report(std::cout);

  Double_ flop, recorded;
  solver.get_flop(flop, recorded);
  std::cout << "\n\nFLOPS= " << (flop/1.0e9)/t
            << ", Modified = " << ((flop-recorded)/1.0e9)/t << "  GFLOPS\n\n\n";

  UHM_ERROR(solver.solve(), ">> Error in Solving");

  Double_ res;
  UHM_ERROR(solver.check(res), ">> Error in Checking");
  std::cout << "\n\nResidual = " << res << "\n\n";

  // This is horribly slow!!!! And even does not reduce residual
  UHM_ERROR(solver.refine(2, res), ">> Error in Refining");
  std::cout << "\n\nResidual = " << res << "\n\n";

  UHM_ERROR(solver.clear(), ">> Error in Clearing");
  UHM_ERROR(solver.flush(), ">> Error in Flushing");

  time_end();

  finalize();

  return true;
}


